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ABSTRACT.  A hardly  known  and  very  important  result  of 
Kreiss  is  proven  explicitly:  Outflow  boundary  extrapolation 

which  complements  stable  dissipative  schemes  for  linear 
hyperbolic  initial  value  problems,  maintains  stability. 

In  view  of  this  result,  the  Lax-Wendroff  and  the  Gottlieb- 
Turkel  schemes  are  applied  to  a test  problem;  as  expected 
from  the  rate-of -convergence  theory  by  Gustafsson,  global  t 
order  of  accuracy  is  preserved  if  outflow  boundary 
computations  employ  extrapolation  of  (local)  accuracy  of  ^ 
the  same  order.  A J 


AMS(MOS)  subject  classification  (1970).  Primary  65M10;  Secondary  65N10 


1.  Introduction 


The  initial  value  problem 


is  well  posed  in  L (Q,°°),  and  requires  no  boundary  conditions  at 
x = 0.  Yet,  it  is  impossible  to  approximate  the  solution  of  (l.l)  by 


a difference  scheme,  which  is  not  right-sided,  without  specifying  boundary 


values  at  some  points  in  a left  neighborhood  of  x = 0 


In  this  paper  we  consider  general  two-sided  dissipative  schemes 


which  are  stable  for  the  pure  Cauchy  problem  for 


main  purpose  is  to  provide  a proof  for  the  following  important  result 
which  was  stated  by  Kreiss  in  1965  [4,  Thm  5],  but  no  detailed  proof 
has  been  published.  We  show  that  if  the  required  boundary  values  are 
defined  by  extrapolation  of  arbitrary  order  of  accuracy,  the  numerical 


algorithm  remains  stable 
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If.  SUPPLEMENTARY  NOTES 


Of  course  one  may  have  instead  of  (l.l),  the  equation  u.  = au  with 
a < 0,  which  defines  a well  posed  initial  value  problem  in  the  quarter 
plane  (x  < 0,  t > 0},  rather  than  in  {x  > 0,  t > 0}.  However,  by 
employing  the  transformation  x -♦  -x  it  is  clear  that  this  problem  goes 
over  to  the  previous  one.  Consequently,  we  would  find  that  the  process  of 
extrapolating  to  grid  points  at  some  right  neighborhood  of  x = 0,  is 
stable. 

To  summarize,  our  aim  is  to  show  that  by  using  a stable  two-sided 
dissipative  scheme  together  with  an  outflow  extrapolation,  to  approximate  a 
well  posed  initial  value  problem  in  the  proper  quarter  plane,  overall 
stability  is  maintained.  Again,  since  the  cases  a > 0 and  a < 0 are 
analogous,  it  suffices  to  prove  stability  for  difference  approximations  of 
(l.l),  and  the  proof  is  given  in  Section  2. 

The  tool  by  which  we  carry  out  the  analysis  is  Kreiss’  stability  theory 
for  dissipative  finite  difference  approximations  of  mixed  initial  boundary 
value  problems.  This  theory  is  given  in  [5],  and  we  assume  that  the 
reader  is  familiar  with  this  work. 

In  Section  3 we  present  numerical  evidence  to  support  the  theoretical 
results.  We  use  two  dissipative  approximations:  The  well  known  centered 

3-point  Lax-Wendroff  scheme,  [7],  and  a centered  5-point  scheme  by  Gottlieb 
and  Turkel,  [2],  In  particular,  our  computations  verify  that  by  using  extra- 
polation of  local  order  of  accuracy  which  equals  the  global  order  of  accuracy 
of  the  difference  scheme,  the  global  accuracy  is  preserved.  The  important 
question  of  convergence  rate  for  mixed  initial  boundary  value  problems  is 
discussed  by  Gustafsson,  [3]. 

The  computations  reported  in  this  work  were  done  on  the  IBM  360 


machine,  at  the  Campus  Computing  Network  of  the  University  of  California, 
Los  Angeles. 


r 


* 
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2.  Stability  analysis 

In  order  to  solve  the  inivial  value  problem  (1.1)  by  a finite  dif- 
ference scheme  let  us  introduce  a mesh-size  Ax  > 0,  At  > 0,  such  that 
X = At/Ax  = constant,  and  use  the  standard  notation  xy  = v Ax,  vy(t) 
v(xy,t).  Now  consider  a dissipative  consistent  approximation  to  (l.l) 
of  the  form 


(2.1) 


vy(t  + At)  = Qvy(t)  , v = 1,2,-..  , 


where 


(2.2) 

and  initial  values  sire  given  by 


if  i 

Q = £ a EJ  , 

j=-r  J 


i-i 


Evv  = vv+l  > 


fi 


(2.3) 


vv(0)  = fv  , v = l,2,...  . 


Here  the  fixed  coefficients  a.  depend  on  a and  X,  such  that  a , a 

j -r  P 

do  not  vanish. 

The  assumption  of  dissipativity  (in  the  sense  of  Kreiss)  means  that 
there  exist  a constant  6 > 0 and  a positive  integer  u,  so  that  the 
amplification  factor 


(2.4) 

of  the  difference  scheme,  satisfies 


§(t)  = £ a.e1^  , -7T  < £ < tt  , 

j=-r  J 


(2.5) 


l$(l)l  < 1 - «lll  , N/UI  < tt  . 


Condition  (2.5)  guarantees,  of  course,  the  (strong)  stability  of  the  approxi- 
mation, should  it  be  applied  to  the  pure  initial  value  problem  for  -“  < x < “. 


1 


Our  final  assumption,  as  indicated  in  the  introduction,  is  that  the 


r > 0,  p > 0.  In  fact,  having  dissipativity: 


our  scheme  must  be  two-sided  if  we  simply  require  that  the  difference 


operator  Q will  be  consistent  with  u = a u for  an  arbitrary  value  of 


a — positive  or  negative.  This  result  is  given  in  Corollary  1 of  [1], 
Since  r > 0,  it  is  evident  that  in  order  to  apply  the  numerical 
approximation  to  (1.1),  we  have  to  specify,  at  each  time  step,  boundary 


tion,  utilizing  the  Lagrange  interpolation  polynomial  of  degree  s - 1 


s > 1,  which  has  accuracy  of  order  s 


In  order  to  comply  with  Kreiss'  formulation  in  [5],  we  should  use  the 


inconvenience  which  we  eliminate  as  follows.  Since  there  exists  a unique 


polynomial  of  degree  s - 1 which  coincides  with  a given  function  at  s 
given  points, (2.6)  is  equivalent  to  extrapolating  from  vQ(t), •••,vfl_^(t) 
to  vQ(t),  and  then  from  vQ(t), • • .,v  1(t)  to  v ^(t),  etc.  That  is, 
(2.6)  is  equivalent  to  the  fixed  coefficient  extrapolation  algorithm 


-5- 


and 


Vx)  = (x  " xu+l^x  ■ XH+2>  •**  (X  - Vs)  l **  = °»***»-r+1 


It  is  straightforward  to  verify  that  regardless  of  n 


(2.10) 


"3 


■ G)'-1)d 


+1 


3 = 


So  our  boundary  conditions  are 

s 


(2.ii)  v^(t)  = e (J)(-Dd+1  Vj(t)  ’ >*  = °>--->-r+1  > 


3=1 


(2.12) 


,?„G)(-1)3  vi 


0 , n = (),•••, -r+i  . 


and  we  finally  write  them  in  the  convenient  form 
s 

r 

d=o 

Now  the  approximation  to  (l.l)  is  well  defined,  where  the  error  at  the 

g 

boundary  extrapolation  is  0(Ax  ). 

Next  denote  by  H the  space  of  all  grid  functions  wy,  defined  for 
v > -r,  which  fulfill  the  boundary  condition  (2.12),  and  which  satisfy 


r |w  |2Ax  < « 
v=-r+l 


(2.13) 

Upon  defining  inner  product  and  norm  by 


(2.H+) 


(u,v)  = T uy  vv  Ax  , 
v=-r+l 


w||  = (w,w)  , 


H becomes  a Hilbert  space. 

According  to  these  definitions  we  may  present  our  difference  approxi- 
mations in  the  form 


(2.15) 


v(t  + At)  = Gv(t)  , v(t),v(t  + At)  e H , 


where  G is  a linear  bounded  operator  in  H defined  by  (2.1)  together 


-5- 


Vx)  = (x  - xu+i)(x  - xu+2^  •**  (x  - xu+s)  i » = 0*  •••#-«■!  • 


It  is  straightforward  to  verify  that  regardless  of  ^ 


(2.10) 


= (®)(-l)J+1,  i = l,"*,s 


So  our  boundary  conditions  are 


g 

(2.11)  yt)  = E (5)(-l)J+1v4+j(t)  , n = 0,...,-r+l 
J — 1 

and  we  finally  write  them  in  the  convenient  form 


(2.12) 


*G>1)d  vj = 0 ’ ^ = °»*-»-r+i  • 


Now  the  approximation  to  (l.l)  is  well  defined,  where  the  error  at  the 

g 

boundary  extrapolation  is  0(Ax  ). 

Next  denote  by  H the  space  of  all  grid  functions  wy,  defined  for 
v > -r,  which  fulfill  the  boundary  condition  (2.12),  and  which  satisfy 


(2.13) 


£ |w  | Ax  < “ . 
v=-r+l 


Upon  defining  inner  product  and  norm  by 


(2.14) 


(u,v)  = T uv  v Ax  , ||w|f  = (w,w)  , 


v=-r+l 


H becomes  a Hilbert  space. 

According  to  these  definitions  we  may  present  our  difference  approxi- 
mations in  the  form 


v(t  + At)  = Gv(t),  v(t),v(t  + At)  e H , 


where  G is  a linear  bounded  operator  in  H defined  by  (2.1)  together 


-6- 


with  (2.12). 

We  say  that  the  algorithm  in  (2.15)  is  stable  if  there  exists  a 
constant  K,  such  that 


(2.16)  !!v(t)||  < K||v( 0 )|]  , \/t  = m^t,  emd  v(0)  € H . 


The  above  description  follows  Kreiss'  representation  in  [5],  so  all 
the  results  of  [5]  hold,  and  we  first  rephrase  the  Main  Theorem  of  [5]  as 
follows  : The  finite  difference  approximation  is  stable  if  the  operator  G 

has  no  eigenvalues  z with  |z|  > 1,  z f 1,  and  if  z = 1 is  not  a 
generalized  eigenvalue  of  G. 

The  concept  of  a generalized  eigenvalue  is  discussed  in  Section  1 
of  [5],  and  in  the  remainder  of  this  section  we  shall  show  that  for  our 
problem, the  Main  Theorem  is  satisfied. 

In  order  to  check  whether  a given  z with  |z|  > 1 is  an  eigenvalue 
of  G,  we  consider  the  characteristic  equation  of  the  difference  operator  Q, 


(2.17) 


det 


P 

r a * 

j=-r  3 


3 


By  Lemma  2 of  [5],  equation  (2.17)  with  z ^ 1,  |z|  > 1,  has  r + p 

roots  x 5 r of  them  with  0 < | k | < 1 and  p with  | k | >1.  More- 
over, according  to  the  proof  of  Lemma  7 of  [5],  as  z -•  1 ( |z|  > 1,  z £ l), 
precisely  one  root  k tends  to  1,  and  this  root  approaches  1 from 
inside  the  unit  disc  if  and  only  if  a < 0.  In  our  case,  the  coefficient 
a in  (l.l)  is  positive,  so  no  root  k of  (2.17)  tends  to  1 from  inside 
the  unit  disc.  Hence,  z = 1 is  not  a generalized  eigenvalue  of  G,  and 
it  remains  to  verify  that  z with  |z|  > 1,  z / 1,  is  not  an  ordinary 
eigenvalue  of  G. 


Suppose  Zq  with  |zQ|  > 1,  zQ  / 1,  is  an  eigenvalue  of  G,  with 
a corresponding  eigenvector  g e H.  That  is,  Gg  = zQg , or  more  specifi- 
cally, due  to  the  definition  of  G,  g must  satisfy  the  relations 


(2.18) 


(Q  - zQ)gv  = 0 , v = 1,2, 


and  the  boundary  conditions 
s 


(2.19) 


j=ox;3/ 


Vj  = 0 ' 


H = 0,-1, ♦••,-r+l  . 


Take  the  characteristic  equation  with  z = zQ,  and  let  k^, 
i = 1,  be  all  its  distinct  roots  which  satisfy  | | <1,  each 

with  multiplicity  7^.  We  know  that  there  are  r such  roots,  so 


(2.20) 


T y.=  T 

i=l 


The  most  general  solution  of  the  ordinary  difference  equation  (2.18), 
which  belongs  to  H,  is  known  to  be 

k-1  v 


q \ 


(2.21) 


gv  = r i »i  k v'  Ki 
v i=l  ltel  1,K  1 


V > -r+1 


where  the  r coefficients  a.  are  arbitrary.  We  still  have  to  verify 

1,  K. 

that  the  solution  in  (2.21)  satisfies  the  boundary  conditions  (2.19).  So, 
we  insert  (2.21)  into  (2.I9)  and  after  a simple  rearrangement  we  get 


q 

(2.22)  r T 

i=l  k=l 

= 0 , 


n = 0, 


-,-r+l  . 


In  (2.22)  we  have  a homogeneous  linear  system  of  r equations  for  the  r 

unknowns  0.  . Denoting  the  coefficient  matrix  by  E,  it  can  be  shown, 

i,k. 

by  elementary  column  operations,  that  E reduces  to  a generalized  Vandermande. 
Hence,  det  E is  proportional  to  an  expression  of  the  form 


I 

f 
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n (k.  - k ) 1 j , 

Ll<i<j^q  J 

where  the  are  integers.  Since  k^,  6116  distinct  with 

0 < |k^|  <1,  we  see  that  det  E ^ 0;  thus  the  only  solution  to  (2.22) 

is  the  trivial  one,  namely  cr  = 0.  Consequently,  g of  (2.21) 

K.  V 

vanishes,  which  means  that  we  have  failed  to  construct  a nontrivial  eigen- 
solution  of  (G  - ZQl)g  = 0 that  belongs  to  H.  Hence,  zQ  is  not  an 
eigenvalue  of  G,  and  Kreiss'  Main  Theorem  assures  stability. 

We  have  proven  the  following  result. 

THEOREM  (Kreiss,  1965).  Let  the  initial  value  problem 

u = au  ; a = const. ; x( sig  a)  > 0,  t > 0 ; u(x,0)  = f(x) > 

be  approximated  by  a dissipative,  strongly- stable,  two-sided  scheme,  which 
is  complemented  by  outflow  extrapolation  of  arbitrary  accuracy,  at  the 
boundary  x = 0.  Then,  the  overall  approximation  is  stable. 


(2.23) 


q 

n 


Pi,, 

Ki  “ Kj/ 


We  conclude  this  section  by  demonstrating  the  cases  r = 1 and  r = 2, 
which  cover  all  schemes  of  practical  importance.  In  particular,  r = 1, 
r = 2,  agree,  respectively,  with  two  schemes  which  we  employ  in 
Section  3. 

For  r = 1,  we  have  to  extrapolate  only  at  one  point,  n = 0,  and 
the  characteristic  equation  has  only  one  root  k inside  the  unit  disc.  So, 
an  eigensolution  of  (G  - z)g  =0,  |z|  > 0,  z ^ 1,  must  be  of  the  form 


g = ox 


(2.24) 


v 


v > 0 


-9- 


and  substituting  this  solution  into  the  boundary  condition  (2.19),  we 
obtain  the  single  equation 


(2.25) 


= '3  0 - 0 • 


(2.26) 


det  E = E = r (j)(-l)J  = (1  - k)S  / 0 , 


and  stability  follows. 

When  r = 2,  we  use  at  each  time  step,  v (t), (t),  to  compute 
v0( t) 5 and  then  vQ(t),  • • •, v^Ct),  to  determine  v 1(t).  Equation  (2. 17) 
has  now  two  roots  with  0 < |k^|  < 1,  and  we  distinguish  between  two 

possibilities.  The  first  is  k = = k,  where 


(2.27) 


gv  = VV  + °2  *V  ’ v Z -1 


The  insertion  of  (2.27)  into  (2.19)  yields  the  system  of  two  equations 

(2.28)  r (5)(-l)J[^  + (U  + 3)a2]  . 0 , n = 0,-1  . 

J —0 

Writing  (2.28)  in  the  form  of  Ecr  = 0,  with  or* a (o^)  being  the  transposed 
unknown  vector,  we  find  that 


(2.29) 


det  E = -k”1(1  - k)2s  ^ 0 . 


The  second  possibility  is  k1  f Kg.  Here 


(2.30) 


gv  = CTi K{  + a2K2  ’ V > -1  > 


which  we  substitute  into  (2.19)  to  obtain 


(2.31) 


£ (jm1  (vr v3*)  ■ 0 - 
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The  coefficient  determinant  is 


(2.32) 


det  E = k’1k"1(1  - Kj/d  - k2)2(k1  - «2)  / 0 > 


and  again,  by  the  Main  Theorem,  stability  follows. 


^ _ Numerical  results 


In  this  section  we  consider  the  initial  value  test-problem 


(3.1) 


u,  = u ; x > 0,  t > 0;  u(x, 0)  = sin  2roc  , 

o X — 


whose  analytic  solution  is 


(3.2) 


u(x,t)  = sin  2rr(x  + t)  . 


We  begin  by  writing  down  the  second  order  accurate  Lax-Wendroff 
scheme  (L  - W),  [7],  which  for  the  linear  equation  in  (3»1)>  takes 
the  form 


(3.3) 


_ m m , a 1 \ 

j=fl  Vv(mAt>  ' 


an  = 1 - X , a+1  = !(X-  + X)  , X = At/Ax  . 


We  recall  (e.g.  [8,  Chapter  12]),  that  in  the  above  case  (a  = 1), 


(3d) 


X < 1 


is  a condition  which  assures  dissipativity  and  strong  stability. 

Here,  r = 1;  hence,  to  approximate  (3.1),  we  need  to  specify  only 
one  boundary  value,  v™.  Extrapolating  via  v™, •••,v™. 


we  get,  according 
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to  (2.11),  a boundary  condition 

(3.5)  vS . , 

which  is  of  (local)  accuracy  of  order  s.  By  the  theorem  of  Section  2, 
the  algorithm  defined  in  (3. 3)  together  with  (3.5),  remains  stable,  provided 
(3.4)  is  satisfied. 

For  the  numerical  computations  we  set  an  artificial  boundary  at  x = b, 
b > 0,  where  we  use  values  of  the  analytic  solution.  Nevertheless,  we 
restrict  attention  to  results  in  the  interval  0 < x < 1,  so  we  choose  a 
large  enough  b,  in  order  to  secure  that  during  the  integration  period, 

0 < t < 1,  there  will  be  no  interaction  between  the  boundary  at  x = b, 
and  the  numerical  solution  at  0 < x < 1.  In  other  words,  errors  due  to 
the  right  boundary,  which  propagate  inward,  never  reach  the  region  0 < x < 1. 


web 

Ax 

s 

■Ham 

r mm 

.05 

2 

20 

2.75  - 2 

40 

5.57-2 

.025 

2 

40 

6.93-3 

80 

1.39-2 

.0125 

2 

80 

1.72-3 

160 

3.46-3 

.05 

B 

20 

OJ 

1 

VO 

• 

Lf\ 

40 

7.03-  2 

.025 

B 

40 

1.98-  2 

80 

2.23  -2 

.0125 

H 

80 

7.00-3 

160 

7.44  - 3 

Table  1:  L-W  results  with  X = At/Ax  = l/2. 

m = t/At  is  the  number  of  time  steps; 

a — n presents  a • 10-n  . 


A 
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The  quantity  l|e||.  in  Table  1 is  the  H-norm  of  the  error,  restricted 

V 0,  ^ ) 

to  0 < x < 1,  i.e.,  in  our  case  (r  = 1), 

(3.6)  !lel!(o,i) * l,e(t)l^0, i)  = SQ[vv(t)  - u(xv^)32ax^  J = l/A*  • 

The  integer  s indicates  the  order  of  extrapolation:  s = 1 and  s = 2 

mean  constant  and  linear  extrapolation,  respectively.  We  realize  that 
all  the  results  are  stable.  As  expected,  linear  extrapolation  maintains 
the  overall  second  order  accuracy  of  the  L-W  scheme,  while  constant 
extrapolation  reduces  the  total  accuracy. 

Our  second  check  relates  to  a centered  5-point  scheme  suggested  by 
Gottlieb  and  Turbel  (G-T),  [4],  We  consider  the  family  of  schemes  in 
(2.4)  of  [2],  set  its  parameters  to  be  a = l/2,  a = 1,  and  linearize. 

The  approximation  we  get  for  equation  (3.1)  is 


(3.7) 


nH-1 


m 


‘O*1 


V * x(x  * §)>  " ■ id 4 §) • x ■ s 

It  was  shown  in  [2],  that  (3.8)  is  stable  if  and  only  if 


(3.8) 


X < s/2/2  , 


and  if  we  somewhat  sharpen  this  condition  and  require 


(3.9) 


x < */£/ 2 , 


we  have  dissipativity  as  well. 

Unlike  the  L-W  case,  which  is  of  seccnd  order  of  accuracy  both  in 
time  and  space,  the  G-T  approximation  is  of  second  order  of  accuracy 
in  time  and  fourth  order  accuracy  in  space.  By  this  we  mean  that  the 
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truncation  error  e satisfies 

(3.10)  e=  At[0(AtL  ) + 0(At»Ax2)  + (KAx1*)]  < At[0(At2)  + . 

Such  schemes  — see  also  the  Kreiss-Oliger  approximation,  [5]  — have 
advantages  when  dealing  with  problems  whose  solutions  have  strong  space 
variations  but  vary  slowly  in  time.  In  particular,  this  concept  fits 
problems  which  approach  a steady  state.  These  ideas  were  discussed  in  [6], 

In  the  G-T  case  r = 2,  so  we  need  to  specify  vQ(t),  v x(t).  Again, 
by  (2.11)  we  have 

(3.11)  'tf’ji.  (aS)(-1)3+1  V*  ’ 

where  it  is  understood  that  the  evaluation  of  v®  precedes  v®^. 


Ax 

X 

s 

m 

■.5 

l'el'(o,i) 

m 

.1. 

Ile,l(0,i) 

.05 

• 5 

4 

20 

9.42-3 

40 

1.85  - 2 

.025 

.25 

4 

so 

5.75-4 

160 

1.16-  3 

• 05 

• 5 

3 

20 

2.14-2 

40 

2.42  - 2 

.025 

.25 

3 

80 

1.78-  3 

160 

1.92  -3 

Table  2 : G-T  results . 


Table  2 shows  that  the  results  are  stable.  Here,  in  analogy  with 
the  previous  case, 

? J 

(3.13)  IM1(o,i)  - s tvv(t)  " » J = . 

Since  G-T  is  of  fourth  order  accuracy  in  space,  it  is  expected  that 
utilizing  a cubic  extrapolation  (s  = 4),  will  maintain  the  4-space 
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accuracy. As  shown,  results  for  Ax  = .05  and  At  = .025  are  compared 
with  those  for  Ax/2  and  At/it.  Indeed,  the  error  is  reduced  by  a factor 
of  16.  This  ratio  is  destroyed  if  we  use  quadratic  extrapolation. 

As  a final  point  of  rerenence  to  the  figures  in  both  tables  we  mention 


that  the  norm  of  the  numerical  solution  was  for  all  times  0 < t < 1, 
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